function [F, G, obj] = FNMTF(X, rowK, colK)
% Fast Nonnegative Matrix Tri-Factorization 
% Input
%         X: nDim * nSmp
%         rowK
%         colK
% Output
%         G: nSmp * colK
%         F: nDim * rowK
% 
% Optimization problem
%     min ||X - F S G'||^2
%     s.t. F \in {0,1}^{d,m}; G \in {0,1}^{n, c}; 
% 
% Note
% (1). F, G are cluster indicator matrices, Not relaxed other matrices
% (2). S is not constrained to be Non-negative, which means this algorithms
% can be applied when X has negative elements
% (3). The optimization of F, G are decoupled, which makes it is easy to
% optimization. 
% 
% [1]. Fast Nonnegative Matrix Tri-Factorization for Large-Scale Data
% Co-Clustering, Hua Wang, IJCAI, 2011
% 

[nDim, nSmp] = size(X);

% **********************************************************
% Initialize G and F with arbitrary class indicator matrices
% **********************************************************

r = randi(rowK, [1, nDim]);
F = zeros(nDim, rowK);
F(sub2ind([nDim, rowK], (1:nDim)', r')) = 1;

r = randi(colK, [1, nSmp]);
G = zeros(nSmp, colK);
G(sub2ind([nSmp, colK], (1:nSmp)', r')) = 1;

nIter = 200;
obj = zeros(nIter,1);
epsilon = 1e-5;
for iter = 1:nIter
    
    % calculate S by Eq (6) ;
    S = F' * X * G;
    S = diag(1./(sum(F) + eps)) * S * diag(1./(sum(G) + eps));
    
    % calculate G by Eq. (8) ;
    smpCenter = F * S;
    smpIdx = knnsearch(smpCenter', X', 'K', 1);
    G = zeros(nSmp, colK);
    G(sub2ind([nSmp, colK], (1:nSmp)', smpIdx)) = 1;
    
    % calculate F by Eq. (8) ;
    feaCenter = S * G';
    feaIdx = knnsearch(feaCenter, X, 'K', 1);
    F = zeros(nDim, rowK);
    F(sub2ind([nDim, rowK], (1:nDim)', feaIdx)) = 1;
    
    obj(iter) = sum(sum( (X - F * S * G').^2));
    % check convergence
    if iter > 1;
        if abs(obj(iter) - obj(iter-1))/abs(obj(iter)) < epsilon;
            obj(iter+1:end) = [];
            break;
        end
    end
end

if iter == nIter
    warning('Not Converged! More iteration needed');
end

return;
    
